
#!/bin/bash
set -e
tools_path=/mnt/ilustre/app/medical/tools
data_path=/mnt/ilustre/app/medical/.data

source ${tools_path}/.var


if [ -z "$1" ]; then
	echo 'clip_reads.sh <input.bam>'
	exit 1
fi

# using b37.fa as the ref

log=.log
if [ ! -e "$log" ]; then
	:> $log
fi

###--------------- argument may be changed ---------------###



genome_name=b37.fa
echo Using $genome_name as the reference genome
genome_assembly=b37
echo genome assembly is: $genome_assembly
dbsnp_version=138
echo dbsnp version is: $dbsnp_version

data_thread_num=10
cpu_thread_num=4

# java_memory=-Xmx16g
java_memory=16g
echo java memory: $java_memory 2>>$log 1>&2

snpeff_db_version=GRCh37.75 
echo $snpeff_db_version 2>>$log 1>&2

header=\
@HD\tVN:1.4\tGO:none\tSO:coordinate

read_group=\
@RG\\tID:${sample_name}\\tPL:ILLUMINA\\tSM:$sample_name
# For gatk, Each read group must contain the platform (PL) and sample (SM) tags.
# For the platform value, we currently support 454, LS454, Illumina, Solid, ABI_Solid, and CG (all case-insensitive).
# Each read in the file must be associated with exactly one read group.

vcf_path=${data_path}/vcf/gatk/

# there is no space between "=" and "\"
ref_genome=\
${data_path}/ref/b37/$genome_name

echo ref genome is: $ref_genome



###--------------- end argument may be changed ---------------###

###--------------- tools path ---------------###

script_path=${tools_path}/script/
bwa=${tools_path}/bwa-0.7.12/bwa
# bwa=/share/apps/bwa-0.7.5a/bwa
samtools=${tools_path}/samtools-1.2/samtools
picard_path=${tools_path}/picard-tools-1.119/
gatk_path=${tools_path}/GenomeAnalysisTK-1.6-9-g47df7bb/
gatk3_path=${tools_path}
gatk=${gatk3_path}/GenomeAnalysisTK.jar
snpeff_path=${tools_path}/snpEff
snpeff=${snpeff_path}/snpEff.jar
snpsift=${snpeff_path}/SnpSift.jar
bcftools=${tools_path}/bcftools-1.2/bcftools
bgzip=${tools_path}/htslib-1.2.1/bgzip
tabix=${tools_path}/htslib-1.2.1/tabix
fastuniq=${tools_path}/FastUniq/fastuniq
###--------------- end tools path ---------------###



echo 2>>$log 1>&2
echo gatk clipreads 2>>$log 1>&2
java -Xmx$java_memory -jar $gatk \
     -T ClipReads \
     -R $ref_genome \
     -I $1 \
     -o $base.trimmed.bam \
     -XF seqsToClip.fasta \
     -X CCCCC \
     -CT "1-5,11-15" \
     -QT 10 \
	 2>>$log